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Chapter  1 

Edge  Constraints 


We  have  developed  discretizations  that  capture  the  singular  behavior  of  cur¬ 
rents  on  open  edges  in  a  high-order  manner.  These  discretizations  use  the 
analytic  properties  of  currents  along  infinity  long,  open  edges.  Consider  the 
metallization  that  covers  the  half  plane  x  >  0.  The  current  normal  to  the  edge 
at  x  =  0,  ~  x1/2+9,  q  G  [0 . . .),  and  the  current  transverse  to  the  same  edge  is 
proportional  to  x~l^2+q.  The  appropriate  testing  functions  and  quadratures 
are  the  Gauss- Jacobi  polynomials  and  quadratures.  All  experiments  indicate 
that  this  is  a  very  good  method,  where  the  number  of  points/ wavelength  to 
achieve  a  given  accuracy  is  of  the  order  of  what  is  need  for  regular  discretiza¬ 
tions  of  regular  geometries.  The  constraint  that  currents  don’t  go  off  the 
edge  is  implicit  in  the  discretization  of  the  integral  equations. 

The  singular  behavior  of  current  along  creases  is  much  more  complex. 
However,  the  current  along  the  edge  is  infinite,  with  leading  singularity  pro¬ 
portional  to  1  where  </>  is  the  angle  of  the  bend,  while  the  current  across 
the  crease  is  finite.  It  is  only  very  special  geometries  with  special  excita¬ 
tions  that  the  error  in  the  current  across  the  edge  dominates.  The  way  we 
have  chosen  to  capture  the  singular  current  along  the  edge  is  to  increase  the 
density  of  discretization  points  by  tapering  patches.  When  the  large  currents 
along  the  edge  have  been  resolved  to  the  desired  accuracy,  the  currents  across 
the  crease  are  captured  even  more  accurately. 

This  characterization  of  the  current  along  geometric  singularities  is  con¬ 
firmed  with  experiments  using  the  automatic  patch  refinement  code  devel¬ 
oped  as  part  of  the  contract.  To  be  specific,  it  is  seen  that  along  the  edges 
of  a  cube,  the  patches  must  be  tapered  towards  the  edge  to  resolve  the  cur¬ 
rent  along  the  edge,  and  this  current  is  orders  of  magnitude  larger  than  the 
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current  flowing  across  the  edge. 

The  singular  discretizations  have  been  applied  to  rectangles,  sharp  elbows, 
tee  structures,  cross  structures  and  PECs  separated  by  gaps  of  width  A/50, 
all  with  good  results.  Care  must  be  taken  when  applying  the  method  to  sharp 
corners  to  ensure  the  singular  currents  are  represented  with  the  discretization. 

Below  are  tables  from  the  simplest  geometry  studied  -  a  3A  x  4A  PEC 
rectangle  with  zero  thickness.  The  solution  used  for  comparison  is  the  12th 
order  singular  discretization.  Using  the  singular  discretization  with  10  points 
per  wavelength,  the  error  is  1.3  •  10-3;  this  compares  well  with  the  accuracy 
obtained  with  a  regular  discretization  of  a  regular  geometry  using  the  EFIE 
and  10  points  per  wavelength.  The  second  table  demonstrates  the  perfor¬ 
mance  of  a  first  order  discretization  of  the  same  rectangle. 


Order 

Points 

Points/ A 

Error 

2 

48 

2 

4.4 

4 

192 

4 

8.6  ■  HT1 

6 

432 

6 

6.1  •  10-2 

8 

768 

8 

5.6  •  10~3 

10 

1200 

10 

1.3  ■  10-3 

Table  1.1:  RMS  error  of  Gauss- Jacobi  method  on  3 A  x  4A  rectangle. 


Points 

Points/ A 

Error 

1200 

10 

9.9  •  KT1 

4800 

20 

3.1  •  HT1 

7500 

25 

2.5  •  HT1 

Table  1.2:  RMS  error  of  first  order  method  with  uniform  mesh  on  3A  x  4A 
rectangle. 
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Chapter  2 

Wave  Quadratures 


2.1  Theory 

The  idea  is  that  in  integrating  oscillatory  functions,  one  can  do  better  than 
Gaussian  quadrature.  The  prescription  is  to  use  formula  that  give  accurate 
results  for  some  spectrum  of  waves,  rather  than  a  finite  set  of  polynomials. 
We  first  illustrate  how  this  can  be  done  on  an  interval. 


2.1.1  Quadrature  on  a  Line 

The  prescription  is  to  pick  weights  (or  “charges”)  wn  and  abscissae  (or  “po¬ 
sitions”)  xn  so  as  to  minimize 


E(x,w,a )  =  f  dfiu(fi)  iune<A1Xn  -  i  [ 

J  n  A  J- 1 


—  J 


EWn^Xn  ~ 


sin(^) 


dxe^x 

2 


=  W{a)  +  Y^wnV(xn-a)  +  wnwn'U (xn  -  xnr,o)  ,(2.1) 

n  n,n' 

where  a  is  a  parameter  of  the  spectrum  w(/x),  and 


W(a)  =  J  dilute)  =  j  y  *  dx  V  (x\  a)  (2.2) 

V{x\a)  =  —2  J  dji  a)  cos  fix— — 
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U(x  —  x'\o ) 


(2.3) 
'  ,  (2-4) 


=  -  f  dx'U{x-x']a) 

=  J  dfiu{ii\<j)e''^x~x''> . 

We  see  the  problem  is  equivalent  to  finding  a  minimum  energy  configuration 
of  particles  with  charges  wn  in  an  external  potential  V (xn)  and  interacting 
with  a  potential  U(xn  —  xn>).  Another  possible  choice  is  a  top-hat  spectrum, 

uj(h)  =  G(a2  —  /z2) ,  (2.5) 

however,  we  chose  to  simplify  the  integrals  by  studying  the  case  of  a  Gaussian 
spectrum. 


Gaussian  Spectrum 

If  we  choose 

then 

W 

V(x) 
U(x  -  x') 


u(p)  = 


\[ko 


-e 


lie 


4-  y/nerf  a 


-V5F/erfi(l  +  x)+erf2(l-x) 


e-V(x-i')2 


(2.6) 


(2.7) 

(2.8) 

(2.9) 


Small  a  limit  If  one  expands  E  in  a  power  series  in  cr2,  one  finds  that  the 
minimization  of  the  successive  term  in  the  series  enforce  the  conditions 

l  ri 

YsWmXm  =  X  /  dxxn ,  (2.10) 

m  L  J  - 1 

so  that  the  small  a  limit  reproduces  Gauss- Legendre  quadrature,  as  should 
be  expected. 


Large  a  limit  In  this  limit,  the  one-particle  potential  V  becomes  very  flat 
except  near  the  endpoints.  The  two-particle  potential  repels  the  points  to 
equal  spacing,  so  that  the  abscissae  become  equally  spaced. 
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Fast  evaluation  of  V ( x )  Although  we  have  a  closed  form  for  V  (x),  in  more 
complicated  cases  (e.g.  triangles),  this  is  not  the  case.  In  such  cases,  a  fast 
evaluation  method  may  be  used  to  make  the  derivation  of  quadrature  rules 
more  efficient.  The  one-particle  potential  is  proportional  to  a  convolution  of 
of  the  two-particle  potential  with  the  domain  of  integration: 


—V(x-,a)  =  J  dx' e  ( 
2  ra/2 

=  ~  d£e~ 
O  J-a/2 


(ax/2-O2 


(2.11) 

(2.12) 


By  using  the  generating  function  for  the  Hermite  polynomials  [AW95],  we  can 
write  this  as  a  power  series  about  the  point  2x  =  £0 /a: 


where 


O  oo 

-V(x;a)  =  -  J2  Cm  fa,  a)  (ox/2-^0)7 

°  m—0 


1  Z^/2 

Cm((o,cr)  =  — 7  /  d£hm{Z-$ o)  (2.14) 

771.  J  ■—& /2 

=  A  [hm-i  {~a/2  -  e0)  -  {a/2  -  &)]  ,  (2.15) 


(2.13) 


=  ~erf({)  (2.16) 

hm(0  =  e-?Hm(Z)  ;  m>0  (2.17) 

and  the  Hm  are  Hermite  polynomials.  In  more  general  cases  the  closed  form 
for  the  coefficients  will  not  be  available.  In  any  case,  they  will  be  evaluated 
for  a  dense  enough  set  of  £0  (these  are  FMM  group  centers)  that  the  power 
series  will  converge  reasonably  rapidly.  In  particular  if  we  let 

.  .  <7  2n  —  1  —  N 


&->{»=  2  N  ;  n  -  1, ....  N , 

(2-18) 

with 

(2.19) 

then 

ox  1 

~2  ~^n  ~  2  ’ 

(2.20) 
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if  we  choose  n  =  \N(x  +  l)/2]  (or  1,  if  x  =  —1).  (The  smallest  integer  not 
less  than  x  is  denoted  by  fx"|.  The  power  series  may  be  truncated  at  a  value 
of  m  equal  to  the  number  of  digits  of  precision  desired.  For  a  1,  direct 
evaluation  of  the  integrals  for  the  coefficients  Cm^n)  becomes  expensive.  To 
deal  with  this,  we  again  use  the  Hermite  generating  function  to  write 


1  fa/2 

crn(£ 0,0-)  =  — r  /  d£ hm  (£  —  £0) 

772!  J~af2 

(2.21) 

(2.22) 

We  then 

use  the  identity 

-&)  =  £  ~hm+l  (£„  -  $>) 

1=0  L 

(2.23) 

to  get 

j  N  oo 

(£o>  &)  =  ~~7  ^2  dnlhm+l  (£n  ~~  £o)  j 

n= 1  1=0 

(2.24) 

where 

dn,  =  Jf  / 

(2.25) 

2  /  a  y+1 

(/  +  !)!  V2ivJ  ' 

(2.26) 

The  last  equality  is  one  for  which  we  have  hope  of  analytic  evaluation 

in  more 

general  cases.  Note  that  the  dni  are  independent  of  £0>  which  is  what  makes 
this  method  “fast” .  The  sum  over  l  may  also  be  truncated  at  the  number  of 
digits  of  precisions  desired. 


Numerical  Minimization 


For  now,  work  directly  with  the  xs  and  ws.  Depending  on  the  minimization 
method  used,  we  may  require  up  to  second  derivatives  of  E  with  respect  to 
the  parameters.  These  are  given  by 


dE 

dwn 


V (xn)  +  2J2  v>n'U  (: xn  -  xn>) 

n' 


(2.27) 


7 


dxn 

^7 iV  (xn)  +  2wny^2)WntU  (xn  3?n') 

(2.28) 

d2E 

71 

—  2 U  ( Xn  ^n') 

(2.29) 

dwndwn’ 

d2E 

^nn'^  (^n)  2'W^tJJ  •^n/) 

dwndxn' 

”f‘25nTj/  ^  ^  WnnU  {xn  Xnf ') 

„ // 

(2.30) 

d2E 

n 

(^n)  ”1“  2 SnnfWn  ^  ^  Wn/fU  (xn  3?n") 

nn 

dxndxn> 

2wnWn>U  {xn  3?n')  • 

(2.31) 

Since  the  energy  is  quadratic  in  the  weights,  they  can  be  eliminated  using 
dE/dwn  =  0.  Letting  E(x)  denote  the  minimum  of  E(x,  w),  with  respect  to 
the  ws,  we  have 


E  = 

TLTif 

(2.32) 

S5 

3h  tqi 

il 

“  o  ^  Unn'Vn'  +  ~pVn  ^2  U U WU *tVnt  , 
n'  Z  Wn’ 

(2.33) 

where 

Vn  =  V(XU) 

(2.34) 

e 

III 

(2.35) 

Unn>  =  U(xn-X'n) 

(2.36) 

Kn>  =  U'(X  n-x’n) 

(2.37) 

EumurJ  =  snn,. 

1 

(2.38) 

Note  that  U'  is  antisymmetric. 

2.1.2  Quadrature  on  Symmetric  Regions 

We  classify  sets  of  group-equivalent  points  in  d  dimensions  and 

parame- 

terize  their  coordinates,  x  —  {aq, . . . ,  Xi, . . . ,  Xd}.  There  are  many  differ¬ 
ent  types  of  indices  in  the  notation  and  it  takes  some  care  to  keep  them 
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straight.  We  take  A  sets  of  abscissae  labeled  by  a,  parameterized  by  {Aa}  = 
Aai,  •  •  ■ ,  Kk,  ■  ■  ■ ,  AjMa  (fie  may  vanish),  having  multiplicity  Ma,  each  asso¬ 
ciated  with  a  weight  wa.  (In  the  cases  we  have  done  so  far,  d  =  2,  so 
0<  fii<  2.)  Thus  for  each  Aa,  there  are  functions  xH(Xa)  for  1  =  1,...,  Ma 
and  i  =  l,...,d.  The  energy  function  to  be  minimized  is  then 

A 

E  =  W(a)  +  £u,aMQI/[Xl(  Aa);a] 

Oc=  1 

A  A 

+  E  E  E  W0U  [xi  (A*)  -  X;  (A^) ;  a]  .  (2.39) 

a=l  0=1  1=1 

[If  we  need  to,  we  could  assume  a  rotationally  symmetric  spectrum,  so  that 
y(x)  =  V(x )  and  U(x)  =  U{x).\  Minimizing  with  respect  to  wa)  we  have 

A  Mp 

0  =  MaV  [x2  (Aa) ;  a]  +  Ma  ^  ^  wpU  [xi  (Xa)  —  x*  (A/?) ;  a] 

(3=1  1=1 

a  mq 

+  E  W0M0  E  U  [X1  (A/?)  -  x/  M  ;  cr]  (2.40) 

P=i  i=i 

Because  U  is  an  even  function,  this  is  a  symmetric  system  of  equations.  In 
fact,  the  two  terms  in  the  U  are  equal  because  the  symmetry  of  the  integration 
is  also  a  symmetry  of  U.  This  can  be  simplified  if  we  define 

M0 

U(Xa,\p;a)  =  Ma  E  U  [xi  (AQ)  —  x;  (A^) ;  a]  =  UaP  =  UPa  (2.41) 

/=i 

V(Xa-,a)  =  MaV  [x2  (Aa) ;  a]  --=  Va .  (2.42) 

Then 

A  _  A 

E  =  W  +  E  waVa  +  E  waUaPwp .  (2.43) 

Q=1  a, 0=1 

As  before,  we  can  either  treat  the  weights  as  independent  variables  or 
eliminate  them. 


Independent  weights 

If  we  consider  the  weights  and  abscissae  as  independent,  useful  partial  deriva¬ 
tives  are 


dE  -  ^  _ 

~  ^  +  2  E  U<*0W0 
°Wa  0=1 


(2.44) 
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where 


dE 

dXaK 

cPE 

dwadi up 
<92£ 

dwad \pK 


&E_ 

d\aK  dXffi, 


A 


+  2wa  E  ^L/3^/3 
/?=1 

(2.45) 

2Ua0 

(2.46) 

A 

fiapVpK  —  2wpU'aKp  +  25a/?  ^  UcikP 

0'=1 

(2.47) 

v4 

8aPwaw0V”nu  +  26a0wa  E 

P’=\ 

-2 ‘WaWpU™,, 

(2.48) 

V' 

v  QK 


77  7 

uclkP 


dVQ 


=  M*E 


dxu(XQ) 


dXQK  ' fr{  dXc 

dU  (XaM 
dX/y* 


Vi  [xi  (Aa) ;  cr] 


M0  d  (\  \ 

E  E  ~5T  a  Vi  [xi  (Aa)  -  X/  (A/j)] 

1  =  1  i  =  1  OACtK 


V" 

r  qku 


dV' 

^  v  CtK 

d\n,v 


=  M,£ 


‘  d2x 


1  i 


-^[xj  (Aa)] 


Tj{2a) 

^awy 


dXaKdXc 

+M„  £  (*«)] 

ij= i  O^otK  O^av 

M~r  d  (  <92Tl. 

Ma  E  E  [X,  (Aa)  -  X/  (A7)] 

l-l  i—i  \  UAaKUAai/ 

,  dxu  (Ag)  ^  ftrij  (Ag) 


n-(2f>) 

^a/c/3^ 


K[x] 


5Aqk 


5Ar 


i=i 

M.EE%^^<j[x,W-x,M]  (2.53) 

(=1  tj=l  OA@u 

dV{x ) 


(2.49) 


(2.50) 


(2.51) 


£4,[xi(Aa)-x/(A7)]J  (2.52) 


OXi 


(2.54) 
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$  ff 
§ 

III 

’x’ 

(2.55) 

^•[x]  = 

(2.56) 

uiA*\  = 

(2.57) 

Dependent  weights 

If  we  eliminate  the  weights  by  minimizing,  we  have 

w.({A})  =  -ii;e-s1({A})V's.  (2.58) 

Z  0=1 

It  is  important  to  remember  that  while  Uap  is  a  function  of  AQ  and  Xp  only, 
U~q  is  a  function  of  all  X.  We  then  have 

*  a, 0=1 
1  A 

w+oJ2  ™eya  (2.59) 

1  Q=1 

-\vLY,u~lpV0  +  \  E  vau-^u'1KQ,u-}pVp 

P  Z  a,a',/3=l 

A 

VL™a  +  2u)a  E  U'aK0™0  ,  (2.60) 

0=1 

2.1.3  Quadrature  on  a  Triangle 

The  quadrature  rule  we  want  approximates  integrals  over  a  triangle,  with 
vertices  chosen  conventionally  to  be 

vi  =  (1,0)  (2.61) 

v2  =  (—1/2,  \/3/2)  (2.62) 

v2  =  (—1/2,  — V3/2) .  (2.63) 

Then  the  integral  we  are  trying  to  approximate  is 

—  f  d2xe-,q  x  =  — —  e*^Sl-S2l/3  — i_e*(42-«3)/3  _ — e*(s3-si)/3  (2  64) 

A  J A  SiS2  s2s3  S3Si  ’  v  ‘  ’ 


E({X})  = 

dE_ 

dXaK 
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where 


A 

III 

CO 

4^  ^ 

to 

(2.65) 

Sl 

=  q-(v2-v3) 

(2.66) 

s  2 

=  q-^s-Vi) 

(2.67) 

=  q  •  (vx  -  v2)  . 

(2.68) 

Gaussian  Spectrum 

The  one  body  potential  is 

V(x;a)  : 

=  -iid2x'u(x-x') 

(2.69) 

where  the  two  body  potential 

is 

U (x;  a)  =  e  °*  x* . 

(2.70) 

Evaluation  of  V(x)  We  have  not  be  able  to  find  a  closed  form  for  the 
single-particle  potential,  but  can  reduce  it  to  one  numerical  integration  (ac¬ 
tually  six  integrations).  This  may  be  sufficiently  fast  that  one  won’t  have  to 
take  the  trouble  to  do  it  really  fast.  The  first  step  is  to  change  the  integral 
over  the  triangle  A  to  an  integral  over  its  boundary  <9 A.  Note  that  we  can 
write  the  two  body  potential  as 

U(x)  =  V2Y(x) 

(2.71) 

where 

Y(x;a)  =  ^ 
a2 

H")+i°g(^)] 

(2.72) 

Then  using  Gauss’  theorem, 

V(x-a)  =  - 

f/sA«ix'n.W(x-x'). 

(2.73) 

Consider  for  the  moment  the  contribution  of  the  integral  over  one  side  (side 
three,  say,  connecting  vj  and  v2)  of  the  triangle  to  this  potential.  With  a 
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suitable  definition  of  h,  l\,  and  l2,  this  is 


VSW  = 


2  2 hr1'  1  -  exp  [-^(/i2 +  x2)] 

~A~^Jo  dX  h2 ~+72 

2  2h  ra2/4  rh  w,2,  2s 
~~A~a2  Jo  dXl  dxe  *  +  [/l^  (2‘76) 


(2.74) 
+  ~ *  h]  (2.75) 


2  yfth  f^!4  e' 

uA 


Z  y/TTfl  f° 

A  a2  Jo 


erf  (VXli) 


Vx 


+  [h  - >  ^2]  (2.77) 


2  .  crli  ah 

— 7  -w  erf —  erf  — 


v/tFT,  r^2/4  ^  e~A;i  erf(y/A/i)1 

Jo 


a2  jo 


dX 


Vx 


J 


2  2  /cr/i  crZj 


+  Ri  — *•  Z2]  (2.78) 

(2.79) 


where 


71 1 <K- 
^1: 


1  -  e-^2 


t/2  +  C2 

e-W 
dX  - 


er: 


f  \/A£ 


VX 


=-A«2 


erf  \/A?7 


VX 


(2.80) 

(2.81) 

(2.82) 


We  expect  the  second  expression  to  be  more  useful  for  77  >  £  and  the  last  for 

£  >  v- 

The  suitable  definitions  are 


h  = 
h  = 

/i  = 


(x  —  Vl)  •  (v2  —  Vl) 
|vi  -  v2| 

(x  -  V2)  •  (V!  -  V2) 
|vi  -  v2| 

\/(x  -  Vj)2  -  T\  . 


=  ^[(x-vi)-(v2-Vl)] 
=  |vx  -  V2|  -  Zj 


(2.83) 

(2.84) 

(2.85) 
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2.2  Numerical  Results 


In  addition  to  the  Gauss-Legendre  rule,  we  have  implemented  two  new 
quadrature  rules  in  the  FastScat  program.  The  first  is  optimized  to  approxi¬ 
mate  the  integral  of  functions  containing  a  fiat  spectrum  of  wavenumbers  on 
a  triangular  patch,  up  to  a  cutoff.  The  second  is  optimized  for  a  Gaussian 
spectrum  of  wavenumbers.  These  quadrature  rules  are  hereafter  called  “Flat 
Wave”  quadrature  and  “Gaussian  Wave”  quadrature,  respectively. 

In  a  direct  comparison  of  these  rules’  ability  to  integrate  sinusoidal  func¬ 
tions  of  different  wavenumbers  on  a  triangular  patch  with  size  of  order  1, 
it  can  be  seen  that  for  low  wavenumbers  (k<ir),  the  Gauss-Legendre  rule 
works  best  regardless  of  how  many  quadrature  points  are  in  the  rule  (Figure 
2.1).  This  is  akin  to  a  heavily  overdiscretized  problem.  For  high  wavenum¬ 
bers  (/c>47r),  none  of  the  rules  converges  for  small  numbers  of  quadrature 
points.  As  the  number  of  quadrature  points  in  the  rule  is  increased,  all  three 
rules  tested  begin  to  converge  (as  measured  by  where  the  slope  becomes  de¬ 
cidedly  negative  in  plots  such  as  Figure  2.2)  at  about  the  same  number  of 
points,  but  both  types  of  wave  quadratures  outpace  the  Gauss-Legendre  rule 
in  reducing  error.  For  intermediate  wavenumbers,  7r<fc<47r,  There  is  no  rule 
with  a  clear  advantage,  and  as  the  number  of  quadrature  points  in  the  rule 
is  changed,  there  are  often  multiple  cross-overs  between  the  amount  of  error 
in  the  different  rules  (Figure  2.3). 

In  tests  scattering  off  of  a  sphere,  both  wave  quadrature  rules  showed 
improvements  in  accuracy  over  Gauss-Legendre  quadrature  under  certain 
circumstances.  For  a  six-patch  sphere  of  radius  A,  the  seven-point  Flat  Wave 
rule  was  about  one  full  digit  more  accurate  than  Gauss-Legendre  across  the 
full  180-degree  bistatic  sweep  (Figure  2.4). 

Any  relative  advantage  of  the  wave  quadrature  rules  over  Gauss-Legendre 
depends  strongly  on  how  converged  the  problem  is  for  the  standard  method. 
In  circumstances  where  the  problem  is  not  well  converged  for  Gauss-Legendre 
quadrature,  using  either  of  the  two  wave  quadrature  rules  does  not  improve 
the  results.  Figure  2.5a  shows  the  result  for  a  radius  2A  sphere  discretized 
in  a  geometrically  similar  way  to  the  1A  sphere  used  above,  and  Figure  2.5b 
shows  the  result  for  another  similar  sphere  of  radius  4A,  by  which  time  there 
is  no  distinction  between  any  of  the  rules  tested,  and  all  are  clearly  not 
convergent. 
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Number  of  Points  in  Rule 


Figure  2.1:  k  —  n/10.  Gauss-Legendre  quadrature  converges  fastest  for  low 
wavenumber.  This  is  the  regime  where  the  Taylor  expansion  of  the  sinusoid 
is  most  valid. 


\ 


Sphere  r=l 


100 

Bistatic  Angle  (Degrees) 


Figure  2.4:  Accuracy  comparison  of  three  quadrature  methods  for  a  sphere 
of  radius  A. 
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Log(0  Relative  Error 


- — -  Gauss-Legendre 
Sphere  r=2X  • — 0  Gaussian  Wave 
■ — •  Rat  Wave 


Sphere  r=4X 


Figure  2.5:  In  regimes  that  are  poorly  converged  for  Gauss-Legendre  quadra¬ 
ture,  wave  quadrature  affords  no  improvement. 


Chapter  3 

Low  Frequency  FMM 


We  developed  a  multilevel  low  frequency  sparse  FMM  operator  for  the  vector 
Helmholtz  kernel  which  has  a  theoretical  apply  time  that  scales  O(N).  It  uses 
a  top-down  automatic  grouping  algorithm,  which  allows  the  groups  to  become 
as  small  as  needed  to  construct  an  efficient  multilevel  operator.  This  operator 
is  implemented  to  run  on  multiprocessor  hardware  with  shared  memory  and 
is  part  of  HRL’s  FastScat  scattering  code. 


points 

Dense 

Time 

LF-FMM 

Time/Scaling 

Max  Rel 
Error 

960 

2.3 

1.2 

5  •  lO"6 

3840 

37 

8.5  /1.44 

1  •  10“5 

15360 

575 

47  /1.23 

1  •  10~5 

61440 

9682 

226  /1.13 

3  •  10"5 

245760 

163432 

993/1.07 

2  •  10~5 

Table  3.1:  Planar  Lattice  of  Spheres.  Apply  time  in  seconds  on  an  Origin 
2000  with  400MHz  R12000/R12010. 

From  these  data  observe  the  asymptotic  approach  to  the  theoretical  linear 
scaling  of  the  multilevel  multipole  operator.  Also  note  the  low  overhead  cost 
of  the  multipole  method  —  the  break  even  point  is  around  one  thousand 
unknowns. 
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points 

Dense 

Time 

LF-FMM 

Time/Scaling 

Max  Rel 
Error 

960 

2.77 

1.14 

8  •  10"6 

3840 

44.4 

9.74/1.55 

7  •  10"6 

15360 

741 

58.5/1.29 

6  •  10~6 

61440 

12043 

452 

3  •  10"5 

Table  3.2:  Planar  Lattice  of  Spheres.  Apply  time/RHS  for  4  RHS;  seconds 
on  an  Origin  200  with  225MHz  R10000/R10010. 


points 

Dense 

Time 

LF-FMM 

Time 

Max  Rel 
Error 

480 

1 

1 

1  •  10“15 

1620 

13 

5 

1  •  10"6 

3840 

71 

32 

3  •  10"8 

7500 

273 

76 

7  •  10~7 

12960 

813 

102 

2  •  10-7 

20580 

2067 

254 

9  •  10“8 

30720 

4936 

576 

5  ■  10"7 

Table  3.3:  Cubic  Lattice.  Apply  times  in  seconds  on  an  Origin  2000  with  a 
400MHz  R12000. 

3.1  Electromagnetic  Multipole  Methods  for 
the  Helmholtz  Equation 

The  field  of  a  current  distribution  can  be  described  in  the  current  free  region 
with  a  multipole  expansion.  Multipole  expansions  are  often  advantageous 
because  complex  current  distributions  can  have  external  fields  accurately 
described  with  a  small  number  of  multipole  coefficients.  When  the  size  of 
the  source  region  is  large  with  respect  to  wavelength  the  external  field  needs 
many  moments  to  describe.  In  this  case,  the  operator  that  translates  the 
multipole  expansions  can  be  diagonalized. 

The  Green  function  for  the  scalar  Helmholtz  equation  is 

I  „ik\x—y\ 

i .  (3.1) 
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For  |x|  >  |y|  this  can  be  expanded  in  terms  of  multipoles  as 
zk 

G(x,  y)  =  ^  E(-)'(2 1  + 1)  W  •  v)  i*(*|y|)  •  (3.2) 

The  dyadic  Green  function  for  the  electromagnetic  equation  is 

1  pik\x-y\ 

G(X,y)  =  ik(  1  +  -VV)^-^.  (3.3) 

Consider  two  non-intersecting  spheres,  each  with  radius  R.  Sphere  X  is 
centered  at  X  and  points  within  the  sphere  are  X  +  x,  |x|  <  R.  Sphere  Y, 
centered  at  Y,  have  points  offset  from  its  origin  by  y. 

For  a  current  distribution  given  by  J(X  +  x)  that  is  zero  outside  sphere 
X,  the  field  E(X  +  x)  outside  sphere  X  can  be  expressed  as 

00  l 

E(X  +  x)  =  J2  £  aim  hi(k\x\)  Ylm(n(x)) ,  (3.4) 

1=0  m=—l 

where  the  vector  coefficients  ct/m  are  determined  by  the  current  distribu¬ 
tion.  The  unit  vector  along  x  is  x.  The  unit  vectors  along  the  cartesian  axes 
are  x ,  y,  z.  The  spherical  angles  are  cos  6  =  x  ■  z  and  cos  </>  =  £■  x,  (  — 
x  —  z  cos  6 .  The  spherical  angles  are  abbreviated  =  (cos  6,<j>)  and  the 
differential  area  of  the  unit  sphere  is  dQ  =  dcp  sin  Odd. 

The  fields  in  a  source  free  region  are  determined  by  the  value  of  the 
fields  on  the  boundary.  Using  the  radiation  condition  at  infinity,  and  the 
orthogonality  of  the  Yim(fl),  the  a/m  are  given  by  the  integration  over  any 
shell  of  radius  A  >  R, 

aim  =  hi\kA)  I  do  Y^(fi)  E(X  +  x(A,  fi))  (3.5) 

The  value  of  A  can  be  arbitrarily  large.  Letting  A  —+  oo, 

aim  =  il+ 'L  f  dn  Y^iQ)  J  d3x  e~ik^>x  [1  +  jfc(ft)  k{ ft)]  •  J(X  +  x) .  (3.6) 
In  sphere  Y,  \y\  <  R,  free  of  currents,  the  field  can  be  expanded  as 

OO  l 

E(Y  +  y )  =  £  £  plm  ji(k\y\)  Ylm(Q(y)) ,  (3.7) 

1=0  m=-l 

The  problem  is  how  to  express  the  in  sphere  Y  in  terms  of  the  alm. 
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3.2  High  Frequency,  kR  >  1 

Write  the  scalar  kernel  between  a  point  x  in  sphere  X  to  a  point  y  in  sphere 
Y  as, 

=  i*£(-)'(2 /  +  1)  P,(S  •  A)  A,(*|A|)  ms\) ,  (3.8) 

where 

5  =  x  -  y,  A  =  X  -  Y,  |A|  >  |*|  (3.9) 

Using  [MF53,  Eq.  11.3.47], 

j,(A:|5|)  Pi{5  •  A)  =  J  dQ  eikf^>5  •  A)  (3.10) 

the  Greens  function  can  be  written 

^l^Al  ‘yT.  . 

^ E*  (2*  +  1)  *«(*|A|)  J  dn  eikfM-5  Pi{f(Q)  ■  A) .  (3.11) 

With  the  truncation  of  the  infinite  sum  to  L  +  1  terms,  the  sum  and 
integration  can  be  exchanged  to  define  a  plane- wave  expansion  of  the  Greens 
function, 

„ifc|5+A[  •£.  - 

—  =  -Jdp.rL(n,A)e‘™>‘,  (3.12) 

where 

L 

Tl(9,  A)  =  J2il(21  +  1)  hi(k\A |)  Pj(r(fi)  ■  A)  (3.13) 

i 

Noticing  that  6  <  2 R,  an  estimate  of  L  comes  from  the  condition  that 

f  dx  Pt(x)  e2ikxR  <  € .  (3.14) 

Prom  this  comes  the  bound,  (k2R)l/l\  <  e;  the  Ith  term  in  the  Taylor 
series  expansion  of  the  exponential  is  less  than  e.  The  empirical  formula 
used  in  FastScat  is 

L  =  k2R  +  -dl^tS  ln(k2R  +  t r) .  (3.15) 

The  choice  of  L  establishes  that  the  order  of  the  quadrature  rule  should  be 
2L,  or  a  Q  =  (L  + 1)2  point  rule  with  abscissae  and  weights,  {9q,  dG?},  given 
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by  a  product  of  a  Gauss-Legendre  rule  in  6  and  a  trapezoidal  rule  in  cp.  The 
Green  function  is  expanded  in  plane  waves  and  is 


p»fc|(5+A|  Ab.  Q 


(3.16) 


3.3  Mid  Frequency 

Why  does  the  high  frequency  kernel  break  as  kR  — >  0?  In  this  limit,  the 
empirical  formula  for  L  is  L  0.7  digits. 

To  get  a  handle  on  how  many  multipole  terms  are  needed,  consider  the 
Laplace  operator  between  the  two  spheres  with  radius  R  centered  at  X  and 
Y  with  a  center  to  center  separation  of  6/2.  Write  the  kernel  l/\5  +  A|.  For 
a  6  colinear  with  A  this  can  be  written 


|<5  +  A|  6/2  1  + 5/6/2  6/2  ^v6/2 


=  —  Y(—y. 

a  r? 


(3.17) 


This  is  valid  because  |5|  <  2/2.  For  this  expansion  to  have  an  error  of  e  for 
the  largest  value  of  8,  L  must  satisfy  the  condition 


or 


(3.18) 

=  -  logs  e 

(3.19) 

which  is,  approximately,  L  =  2  digits. 

Using  this  value  for  L,  the  high  frequency  method  can  still  be  used  if 
&|A|  is  not  too  small.  For  such  values  the  leading  term  in  the  expansion  of 
y/(fc|A|)  is 

V'WA|)  =  <3-20) 

It  is  possible  to  preserve  numerical  significance  in  the  sum  used  to  compute 
A)  if  the  condition 


lfc|A|  > 


l/(2L  —  1)1! 
e  1015 


(3.21) 


holds  true,  for  a  floating  point  precision  of  10-15. 
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The  mid  frequency  method  in  FastScat  uses  the  high  frequency  method 
with  this  modified  value  of  L  and  the  restriction  on  the  minimum  value  of 
fc|A|.  This  method  has  to  be  used  very  carefully  because  it  fails  catastroph¬ 
ically  as  &|A|  goes  to  zero. 


3.4  Low  Frequency 


The  operator  that  shifts  the  apq  of  sphere  X  to  the  fiim  of  sphere  Y  can  be 
determined  by  evaluation  of  the  fields  from  the  sources  in  sphere  X  on  a  shell 
of  non-zero  radius  around  point  Y .  Thus 

film  =  T lm,pq  &pq  5  (3.22) 

r im<pq  =  jr\kB )  J  dft  Y^(n)  hp(kR )  Ypq(Ci) .  (3.23) 

In  this  expression,  R  and  ft  are  the  distance  and  angle  from  point  A  to  a 
point  at  ft  on  the  shell  of  radius  B  centered  at  Y.  The  radius  B  must  be 
such  that  ji(kB)  ^  0.  For  the  low-frequency  domain,  the  choice  B  =  R/2  is 
not  unreasonable.  Using  this  shift  operator  alone  a  single  stage,  O(NlnN), 
method  can  be  developed. 

To  implement  an  0(N )  multilevel  scheme,  the  coefficients  must  be  shifted 
to  and  from  the  centers  of  parent  and  child  groups.  A  child  group  is  contained 
in  the  sphere  centered  at  A  of  radius  R  and  this  child’s  parent  group  is 
centered  at  X  and  has  a  radius  of  2 R.  The  child  group  is  wholly  contained 
within  its  parent  group. 

To  shift  the  apg  centered  at  A  to  a  parent  group  of  radius  2 R  centered 
at  A  near  enough  so  that  |A  —  A|  <  R,  use 

alm  =  hrl(k2R)  J dft  Y^i ft)  hp(kR)  Ypq( ft)  apq  ,  (3.24) 

where  R  and  ft  are  the  distance  and  the  angle  from  point  A  to  a  sphere  of 
radius  2 R  centered  at  A.  To  shift  the  aim  into  the  plane- wave  expansion 
used  in  the  high  frequency  expansion,  use, 

E^ft)  =  32  l-Y  alm  Ylm(ft) .  (3.25) 
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Trickling  down  the  hierarchy,  the  Am  are  expressed  in  terms  of  the  plane 
wave  expansion  as 

Am  =  47 nl  J  dn  VSJQ)  E^Q).  (3.26) 

Within  the  low-frequency  domain,  the  /3pq  are  shifted  down  with 

Pim  =  jr\kR/ 2)  J  dn  Yt*m(n)  jp{kh)  vpq(n)  ppq ,  (3.27) 

where  R~  and  Cl  are  the  distance  and  the  angle  from  the  center  of  the  parent 
group,  X,  to  a  point  at  n  on  the  child’s  shell  of  radius  R/2  centered  at  X. 
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Chapter  4 

Automatic  Patch  Refinement 


We  developed  a  method  to  refine  the  discretization  and  have  validated  the 
method  on  a  variety  of  complex  geometric  objects.  For  each  run  of  the 
scattering  program  the  currents  are  saved  to  a  file.  A  postproccessing  step 
reads  this  file  and  determines  which  patches  have  the  largest  current  error, 
then  it  subdivides  the  patches  with  the  largest  errors  and  generates  a  new 
geometry  file.  The  jet  is  defining  a  suitable  measure  of  current  error,  and 
recognizing  that  a  patch  can  be  split  in  two  directions,  one  of  which  may  be 
adequate  to  resolve  the  current. 

On  each  patch,  two  components  of  the  surface  current  are  known  at  uq, 
the  Q  nodes  of  a  quadrature  rule,  //  e  [1,2]  and  is  the  orthogonal 
tangent  basis.  The  current  at  point  q  is 

jq  =  jl,  q^l(Uq)  +  jl,  q^iUq). 

A  set  of  scalar  functions,  4>k(u),  can  be  fit  to  these  values,  with, 

K 

jp,q  0  (uq)  jp}k  =  *&qkjk- 

k= 0 

If  the  number  of  independent  functions  equals  the  number  of  quadrature 
points,  the  system  is  square.  A  reduced  system  of  functions  can  be  con¬ 
structed  from  the  square  set  by  differentiation  with  respect  to  the  parametric 
directions. 

K 

•A1.?  ^  ^  dy(j)  ( Uq )  j^  f,  =  A qkjk- 

k= 0 
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This  system  of  over-determined  equations  is  singular;  the  coefficients  can 
be  found  in  a  least  squares  sense  using  singular  value  decomposition  of  the 
matrix  A  and  its  pseudo-inverse,  AM.  The  array  of  errors  at  quadrature 
points,  ug,  q  €  [1 . . .  Q]  is  given  by 

e^  =  (l-AA~1)^,  ^  =  1,2.  (4.1) 

A  measure  of  the  error  over  the  entire  patch  is  given  by 

=  M  =  1*2)  (4.2) 

Q 

where  wq  are  the  quadrature’s  weights.  For  convex  scatters,  this  measure  of 
error  corresponds  well  with  the  error  in  the  far  field.  For  other  convoluted 
geometries,  there  is  no  such  simple  measure. 

If  the  error  e ^  is  larger  than  a  user  specified  threshold,  then  the  patch  is 
divided  into  two  halves  (in  parametric  space).  This  leads  to  the  result  that 
a  quadrilateral  patch  could  be  divided  into  2  or  4  sub-patches,  or  not  at  all. 
It  is  worthwhile  to  investigate  the  error  at  each  point  in  the  hope  of  finding 
a  point  more  suitable  than  the  parametric  center  to  divide  the  patch. 

If  there  are  multiple  computed  currents  corresponding  to  different  exci¬ 
tations  (e.g.  resulting  from  an  angle  scan),  then  the  error  is  the  RMS  error 
of  all  the  currents.  As  the  code  is  currently  written,  there  are  at  least  two 
currents  corresponding  to  the  two  polarizations  of  an  incident  plane  wave. 

This  code  has  been  tested  on  a  variety  of  geometries,  with  PEC  and 
penetrable  volumes  The  example  test  case  summarized  in  Table  4.1  is  a 
constellation  of  PEC  objects  (  a  sphere,  a  thin,  flat  plate,  and  a  flattened 
torus,  in  close  proximity  to  each  other  )  with  a  dipole  source  a  small  fraction 
of  a  wavelength  off  the  surface. 
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Pass 

Points 

Error 

i 

1464 

4.17- 10"1 

3 

1626 

1.47- 10-1 

5 

2184 

1.02  •  10-1 

7 

3138 

6.06  •  10~2 

9 

5586 

3.35  •  10"2 

11 

8106 

1.73  •  10“2 

13 

11994 

7.28  •  10-3 

15 

15144 

3.11  •  HT3 

Table  4.1:  Convergence  of  Automatic  patch  refining  algorithm. 
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Chapter  5 


Raytheon  Missile  Systems 
Code  Validation 

A  subcontract  was  awarded  to  RMSC  for  validation  of  the  FastScat  code 
against  range  measurements.  Jack  Kennedy,  an  RMSC  engineer,  performed 
the  work  on  the  subcontract.  The  appendix  of  this  document  contains  a 
complete  report  of  his  investigation. 

The  validation  task  included  construction  and  RCS  range  measurement 
of  three  metal  targets  and  numerical  simulation  of  the  RCS  for  these  targets 
using  three  electromagnetic  modeling  codes. 

The  targets  were  chosen  to  test  the  ability  of  the  code  to  make  accurate 
RCS  predictions  for  targets  containing  geometric  features  commonly  found 
on  real  air  vehicles.  The  first  target,  a  cube  six  inches  on  a  side,  was  chosen 
to  test  the  effectiveness  of  each  code  in  dealing  with  targets  that  contain 
geometric  singularities,  in  this  case,  edges  and  corners.  When  geometric 
singularities  are  present  they  can  interrupt  the  predominant  movement  of  the 
surface  current,  producing  standing  waves  and  strong  radar  reflections.  They 
also  result  in  surface  current  singularities,  which  pose  a  modeling  challenge. 
The  second  target,  a  trihedral  corner  reflector  with  sides  of  4,  6,  and  12 
inches,  was  chosen  to  test  each  code’s  ability  to  predict  the  RCS  of  a  target 
for  which  multiple  bounce  scattering  plays  a  significant  role.  The  RCS  of 
this  target  is  also  affected  by  the  presence  of  numerous  edges  and  corners. 
The  third  target  is  an  8”  long  by  8”  wide  -wing-like  fin.  Its  surface  is  mostly 
smooth  except  for  an  acute  angle  corner  along  the  trailing  edge.  This  target 
was  chosen  to  test  how  well  each  code  computes  the  RCS  of  a  curved  target. 

Each  target  was  take  to  the  Raytheon  Site  A  outdoor  test  range  in  San 
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Diego,  CA  where  monostatic  RCS  measurements  at  multiple  aspects  and 
frequencies  were  obtained.  Details  of  the  measurement  process  can  be  found 
in  the  full  report. 

The  various  measurements  were  compared  against  computed  results  from 
three  electromagnetic  modeling  codes:  FastScat,  FISC,  and  Xpatch.  For  the 
two  targets  comprised  of  flat  surfaces,  the  results  computed  by  FastScat  and 
FISC  matched  the  experimental  results  well.  Both  codes  were  notably  su¬ 
perior  to  Xpatch  in  terms  of  RCS  accuracy,  especially  at  low  frequency.  For 
the  curved  target,  FastScat  and  FISC  were  again  quite  superior  to  Xpatch. 
In  this  case,  however,  FISC  was  generally  less  accurate  than  FastScat,  espe¬ 
cially  for  lower  frequencies.  The  difference  can  be  attributed  to  the  fact  that 
FastScat  internally  represents  a  curved  surface  using  curved  mesh  elements 
whereas  FISC  approximates  a  curved  surface  using  flat  mesh  elements.  For 
this  reason  and  others,  it  was  generally  possible  to  obtain  more  accurate 
results  from  FastScat  in  less  time  than  for  FISC. 

The  overall  conclusion  is  that  FastScat  is  a  significantly  advanced  electro¬ 
magnetic  computation  tool  that  generally  produces  superb  RCS  predictions 
for  metal  targets. 
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Appendix  A 

Fastscat  Validation,  Final 
Report 


This  report  is  supplied  as  a  separate  volume. 


32 


Bibliography 


[AW95]  George  B.  Arfken  and  Hans  J.  Weber.  Mathematical  Methods  for 
Physicists.  Academic  Press,  San  Diego,  fourth  edition,  1995. 

[MF53]  Phillip  M.  Morse  and  Hermann  Feshbach.  Methods  of  Theoretical 
Physics.  McGraw-Hill,  New  York,  1953. 

[WX99]  Stephen  Wandzura  and  Hong  Xiao.  Quadrature  rules  on  triangles 
in  R2 .  Technical  Report  YALEU/DCS/RR-1168,  Yale  University, 
Department  of  Computer  Science,  January  1999. 


33 


